Detection of 6 K gas in Ophiuchus D 
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Abstract 

Cold cores in interstellar molecular clouds represent the 
very first phase in star formation. The physical conditions 
of these objects are studied in order to understand how 
molecular clouds evolve and how stellar masses are deter- 
mined. The purpose of this study is to probe conditions 
in the dense, starless clump Ophichus D (Oph D). The 
ground-state (1 io — * In) rotational transition of ortho- 
H2D + was observed with APEX towards the density peak 
of Oph D. The width of the H2D + line indicates that the 
kinetic temperature in the core is about 6 K. So far, this 
is the most direct evidence of such cold gas in molecular 
clouds. The observed H2D + spectrum can be reproduced 
with a hydrostatic model with the temperature increasing 
from about 6 K in the centre to almost 10 K at the surface. 
The model is unstable against any increase in the exter- 
nal pressure, and the core is likely to form a low-mass 
star. The results suggest that an equilibrium configuration 
is a feasible intermediate stage of star formation even if 
the larger scale structure of the cloud is thought to be de- 
termined by turbulent fragmentation. In comparison with 
the isothermal case, the inward decrease in the tempera- 
ture makes smaller, i.e. less massive, cores susceptible to 
externally triggered collapse. 



Direct measurement of the gas properties, such as 
the kinetic temperature, is possible using spectral lines. 
Spectral line observations of very dense, starless cores 
are, however, hampered by the fact that common tracer 
molecules like CO and CS freeze there onto dust grains 
(Willacy et al. 1998; Caselli et al. 1999). The nitroge- 
nous compounds NH3 and N2H + seem to survive longer 
than most other species (Tafalla et al. 2006). This is for- 
tunate because the rotational temperature of NH3 is con- 
sidered to be a good measure of the kinetic temperature. 
Recently NH3 excitation was used to derive a very low gas 
temperature, ~ 5.5 K, in the centre of the prestellar core 
LI 544 in Taurus (Crapsi et al. 2007). 

At very high densities (n ~ 10 6 - 10 7 cm 4 ) and low 
temperatures (T < 10 K) also NH3 and N2I-L are likely 
to freeze out, and the gallery of useful molecular trac- 
ers becomes very limited (Walmsley et al. 2004). In 
these circumstances the abundances of H3 , and its ener- 
getically favoured deuterated forms, H2D + , D2FL, and 
Dt, are predicted to increase strongly. The asymmetric 
isotopologues H2D + and D2H + have permanent electric 
dipole moments, so can be used as probes of these oth- 
erwise hidden regions (Caselli et al. 2003; Vastel et al. 
2004; van der Tak et al. 2005). 



1 Introduction 

Starless cores of molecular clouds are heated externally 
by the interstellar radiation field (ISRF) and by cosmic 
rays. Theoretical estimates of the attenuation of the ISRF 
by dust suggest that the temperature decreases to about 5- 
6 K in the centres of dense starless cores (Evans et al. 
2001; Zucconi et al. 2001; Stamatellos & Whitworth 
2003). Observational evidence of low temperatures is still 
scarce. 

"This publication is based on data acquired with the Atacama 
Pathfinder Experiment (APEX). APEX is a collaboration between the 
Max-Planck-Institut fiir Radioastronomie, the European Southern Ob- 
servatory, and the Onsala Space Observatory. 



Here we use the singly deuterated trihydrogen ion, 
H2D + , to show that the gas kinetic temperature is about 
6 K in the very dense core of the starless clump Oph D 
(LI 696 A), which lies in the nearby p Ophiuchii molecular 
cloud. The position observed is located near the southern 
tip of the elongated clump. It corresponds to the thermal 
dust emission maximum at A = 850/vm on the SCUBA 
map of Kirk et al. (2005), and the dust absorption maxi- 
mum on the ISOCAM 7/im map of Bacmann et al. (2000). 
This position coincides with the N2D + peak in the survey 
of Crapsi et al. (2005; see their Fig. 10), and represents 
the density and column density maximum of the clump. 
The 3D structure of the whole clump has been recently 
modelled by Steinacker et al. (2005). 
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Figure 1: Left: The H2D + spectrum towards Oph D. Right: ISOCAM 7/jm mid-infrared absorption image of the 
Oph D clump obtained by Bacmann et al. (2000; the image is reproduced from Bergin & Tafalla 2007). The scale is 
indicated by a 1' bar. 



2 Observations and data reduction 

The observations were made with APEX on May 15, 
2006. The telescope is described by Giisten et al. (2006). 
The 372421.364 MHz line of ortfio-H 2 D + was observed 
in the upper sideband of the APEX-2A SIS DSB re- 
ceiver. The HPBW of the antenna is 17" at this frequency. 
The backend was the Max-Planck-Institut filr Radioas- 
tronomie Fourier transform spectrometer (FFTS). The 1 
GHz band of the FFTS was divided into 16384 chan- 
nels resulting in a channel width of 61 kHz which cor- 
responds to ~ 50 ms~' at the observed frequency. The 
effective resolution of the spectrometer corresponds to 80 
ms~'. The observations were performed in the position- 
switching mode. The coordinates of the ON-position are 
R.A.16 h 28 m 28' s :9, Dec.-24°19'09" (J2000). The OFF- 
position was selected -10' south of ON-position in Dec. 
The spectrum shown in Fig. 1 has been obtained by adding 
60 ON-OFF measurements with a duration of 20 seconds 
per phase. The observed position is indicated on the 7yum 
ISOCAM absorption map of Bacmann et al. (2000; repro- 
duced from Bergin & Tafalla 2007). The observing con- 
ditions were good (PVW 0.34 mm, zenith opacity 0.33 at 
372 GHz). Oph D was in the elevation range 60° - 70°, 
and the DSB system temperature was about 210 K during 
these measurements. 

The observed spectra were reduced in a standard way 
using the GILDAS software packag^H- The occasional 
low-frequency ripple in the individual raw spectra was 

'Grenoble Image and Line Data Analysis Software 
package has been developed by IRAM-Grenoble, see 
|http : //www . iram . f r/ IRAMFR/GI LDAS| 



first fit with a sinusoidal baseline, after which possible 
higher frequency ripple was masked out in Fourier space. 
After averaging the 60 individual spectra, a first-order 
baseline was subtracted. 

3 Results 

The detected H 2 D + line has a single, narrow component. 
A Gaussian fit to the spectrum gives the following line pa- 
rameters: = 0.60 ± 0.09 K, v LSR = 3.40 + 0.01 kms" 1 , 
and Av(FWHM) = 0.26 + 0.03 kms" 1 . The molecular 
mass of H 2 D + is only 4 a.m.u. The observed linewidth, if 
assumed to be caused only by thermal broadening, indi- 
cates a kinetic temperature of Jki n = 6.0 + 1.4 K. This is 
an upper limit as the estimate neglects the possible non- 
thermal broadening and the slight instrumental broaden- 
ing, the contribution of which is about 0.01 kms" 1 . 

4 Discussion 

The temperature obtained towards Oph D is similar to 
the estimate for the nucleus of LI 544 mentioned above 
(Crapsi et al. 2007). The present determination is more 
direct than that in LI 544 because it does not involve mod- 
elling. These two results using very different methods 
consolidate evidence of very low temperatures inside star- 
less dense cores. 

The H 2 D + position coincides with the N 2 D + maximum 
of the core with an exceptionally high N 2 D + /N 2 H + col- 
umn density ratio (Crapsi et al. 2005). The enhancement 
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of N2D + is caused by the precursor molecule N2 reacting 
with increasingly abundant H2D + , D2H + or instead of 
the normal isotopologue Hi. Carbon monoxide has been 
observed to be heavily depleted in this core (Bacmann et 
al. 2002). The accumulated observational results, i.e. CO 
depletion, a large degree of deuterium fractionation, the 
detection of H2D + , and a very low kinetic temperature, 
all fit qualitatively into the current picture of chemistry in 
dense starless cores (Walmsley et al. 2004). 

The picture is not complete, however, until we man- 
age to build a physical model for the object to explain 
the observed H2D + profile. The dense central core ob- 
served here apparently has very little turbulence and looks 
roundish in the N2D + maps (Crapsi et al. 2005). It 
may therefore be in near hydrostatic equilibrium. As the 
H2D + line, which probably originates in the core centre, 
indicates a very low temperature, it is reasonable to as- 
sume a radial temperature gradient. We adopt the modi- 
fied Bonnor-Ebert model (Evans et al. 2001; Zucconi et 
al. 2001), i.e., a self-gravitating, hydrostatic core with 
a temperature gradient caused by the attenuation of the 
ISRF. As the condensation is part of a larger cloud, we 
assume, somewhat arbitrarily, that the obscuration by the 
surrounding envelope corresponds to a visual extinction 
of Ay = 10 m . This choice is of little importance because 
the obscuration by the core itself must be much higher to 
produce a temperature close to 6 K in the centre. 

4.1 Core model 

The density and temperature distributions in the core were 
solved in an iterative manner. The density distribution was 
first solved for an isothermal Bonnor-Ebert sphere (Alves 
et al. 2002), and the radial temperature gradient was cal- 
culated using the ISRF and dust opacity models adopted 
from the literature (Black 1994; Ossenkopf & Henning 
1994). This temperature distribution was used to derive a 
new density profile by integrating the equation of hydro- 
static equilibrium. 

We assumed that the gas temperature, 7\i n , is equal 
to the dust temperature, T^ust- This assumption is gen- 
erally believed to be valid at densities relevant to this 
study (> 10 5 cm" 3 , Burke & Hollenbach 1983). Recently, 
Bergin et al. (2006) found evidence of different gas and 
dust temperatures in the centre of the globule B68 with 
n(H2) ~ 3 10 5 cm" 3 . We note, however, that the nearly 10 
times higher density in Oph D is likely to result in a closer 
dust-gas thermal coupling as compared with B68. 

All radiative transfer calculations were made using our 
Monte Carlo code (Juvela 1997). It was found that cen- 
tral densities in excess of 10 6 cm" 3 are needed to make 
the temperature decrease close to 6 K. As a test, we com- 
pared the temperature profile calculation with the results 
of Stamatellos & Whitworth (2003). Identical tempera- 



ture distributions were obtained using the same models 
for the core stucture, ISRF spectrum, and the dust opac- 
ity. 

4.2 H 2 D + excitation 

The structure and rotational spectrum of H2D + are well 
known (Miller et al. 1989). Like molecular hydrogen, 
H2, the molecule has two nuclear spin states, ortho (here- 
after o-H2D + , the H nuclei have parallel spins), and para 
(p-H2D + , opposite spins). While radiative transitions are 
only possible between rotational levels of the same nu- 
clear spin state, collisions with H2 can also result in ortho- 
para conversion (e.g., Gerlich et al. 2002). 

The state-to-state coefficients for the H2D + + H2 sys- 
tem are not available. In the recent study of H2D + pro- 
files towards L1544 van der Tak et al. (2005) used scaled 
radiative rates adopted from Black et al. (1990), where 
the downward collisional coefficients within the ortho and 
para ladders were assumed to be proportional to the line 
strengths of the corresponding radiative transitions, and a 
constant value was used for downward inter-ladder transi- 
tions. The upward rates were evaluated by using the prin- 
ciple of detailed balancing. 

A microcanonical statistical study of the reaction H^ + 
H 2 has been published recently (Park & Light 2007). This 
treatment is being applied to the deuterated isotopologues, 
but the state-to-state reaction probabilities for H2D + + H2 
have not yet been calculated (Hugo et al. 2007). In the 
present paper the collisional coefficients were calculated 
using a modification of the Oka & Epp (2004) model. The 
coefficients are normalized to conform with the available 
chemical data on the nuclear spin-changing reactions. The 
assumptions and parameter values used in the calculation 
are available in the appendix. The approximation corre- 
sponds to the canonical approach discussed in Hugo et al. 
(2007). 

The core model and the collisional coefficients were 
used to calculate the populations of the rotational levels 
of H2D + , and H2D + (lio — > In) line profiles towards the 
core centre. The calculated spectrum corresponds to the 
source brightness temperature distribution convolved with 
the 17" APEX beam. The central density and the outer ra- 
dius of the core and the total (o+p) H2D + abundance were 
varied, and the resulting line profile was compared with 
the observed one. Our best-fit H2D + spectrum and the 
corresponding core model are shown in Fig. 2. The model 
spectrum was converted to the T* K scale by multiplying it 
by the main-beam efficiency, r/ MB = 0.73. 

4.3 Modelling results 

The line intensity can be reproduced by a set of hydro- 
static models with the central density and the central tem- 
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Figure 2: The core model and the calculated H2D + spec- 
trum. Top: Density as a function of radial distance 
from the core centre according to our hydrostatic, non- 
isothermal model (continuous line). For comparison the 
solution for an isothermal sphere with 7\i nj i so =< T^ n >= 
8.3 K and the same peak density is shown as a dashed 
line. Middle: The radial profiles of the kinetic tempera- 
ture, Tkin (assumed to be equal to T^st), and the excitation 
temperature of the line, T ex . The line is thermalized in the 
core centre, but becomes subthermally excited towards the 
edge as the density decreases. Bottom: H2D + (lio -» In) 
spectrum (in the scale) resulting from the Monte Carlo 
radiative transfer calculations applied to the hydrostatic 
model presented in the two upper panels. The observed 
spectrum is shown as a dotted line. 



perature ranging from 2 10 6 cm 4 to 1 10 7 cm 4 , and from 
6.8 K to 6.1 K, respectively. The best-fit model has a cen- 
tral density of 4 10 6 cm 4 . In this model the temperature 
decreases from 9.6 K at the surface to 6.4 K in the cen- 
tre. The line width of the corresponding model spectrum 
is 0.35 kms -1 . This and all model spectra having similar 
peak intensity are broader than the observed one. There 
are two reasons for this: 1) the kinetic temperature in the 
centre of the model is slightly higher than 6 K, and 2) the 
peak optical thickness of the line must be close to unity to 
give sufficiently high intensity. Opacity broadening is re- 
sponsible for about one half of the 'extra' linewidth (~ 0.1 
kms -1 ) of the best fit spectrum with T pea k ~ 1.1. This 
t value results from a total H2 column density of 1 10 23 
cm -2 and a fractional H2D + abundance of 5 10" 10 . The 
core mass is nearly the same, 0.25 M e , for all hydrostatic 
models quoted above. The outer radius was assumed to 
be 3900 AU (apparent radius 30")- This is consistent with 
the published maps (Kirk et al. 2005; Crapsi et al. 2005). 
We assume a distance of 130 pc (Knude & H0g 1998). 

The core can be made cooler by increasing the cen- 
tral density further. This would, however, decrease the 
line intensity because of a larger optical thickness and 
beam-dilution effects as the brightness distribution be- 
comes more centrally peaked. 

We also calculated temperature distributions and H2D + 
spectra for the density profiles derived previously for Oph 
D by Motte et al. (1998) and Bacmann et al. (2000), 
who used millimetre dust continuum emission and mid- 
infrared absorption, respectively. These models assume a 
constant density within a radius of 20" -25" from the core 
centre, and a power law in the outer envelope up to 85". 
For the Motte et al. (1998) model, with a central density of 
9 10 5 cm 4 , the calculated temperature goes down to 6.1 
K in the centre. The Bacmann et al. (2000) model with a 
lower central density (3 10 s cm 4 ) gives a minimum tem- 
perature of 7.4 K. The model of Motte et al. produces sim- 
ilar H2D + spectra to the hydrostatic models quoted above. 
The existence of a large, homogenous, and non-turbulent 
nucleus could possibly be justified by invoking a stabiliz- 
ing magnetic field. 

In view of the moderate S/N of the observed H 2 D + line 
and the uncertainties related to the collisional coefficients, 
the hydrostatic model gives a reasonably good fit to the 
observed spectrum. A high S/N H2D + spectrum towards 
this source, together with rigorously derived collisional 
coefficients will show if the slightly 'too bright and too 
narrow' line remains a problem. Assuming that the excita- 
tion within the ortho and para ladders is thermal, the H2D + 
spectrum constrains the model parameters in the follow- 
ing way. 1) In order to reproduce both the observed inten- 
sity and linewidth, the central temperature must be about 
6 K but not much less. The minimum temperature sets a 
limit to the maximum density through the adopted ISRF 
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and dust opacity models. 2) The H2D + line is preferably 
optically thin, and in any case the peak optical thickness 
cannot be much larger than unity. For high opacities, the 
line peak flattens and the FWHM becomes clearly larger 
than observed. This determines the upper limit of the col- 
umn density and the fractional abundance of the molecule. 
According to the adopted model for ortho-para transitions, 
the o/p ratio of H2D + increases as the temperature de- 
creases, and consequently, in the centre most H2D + is in 
ortho states. This effect facilitates the detection of H2D + 
in cold cores. Furthermore, the accuracy of the modelling 
results is improved by nearly all H2D + molecules being in 
the two states coupled by the observed transition. 

According to Steinacker et al. (2005) the elongated 
structure of the Oph D clump is likely to have formed as 
a result of turbulent fragmentation. They suggest that the 
southern condensation (i.e. the core observed here) may 
be gravitationally bound and collapsing. The present ob- 
servation and modelling conform with the idea of a self- 
gravitating core, but the very narrow line profile suggests 
that the collapse has temporarily halted. 

Compared with an isothermal core at a kinetic temper- 
ature Ikjn iso, a core with inward decreasing temperature 
with the same central density and the average temperature 
< 7\in >= 7kin,iso has a smaller critical radius (Alves et al. 
2002) beyond which the core becomes unstable against an 
increase of the external pressure. While for an isothermal 
sphere the critical dimensionless radius parameter, is 
6.45, the corresponding parameter for our best-fit model 
is 6.20. This value corresponds to 1700 AU (13"), which 
is clearly smaller than the core radius derived in previous 
studies. Even though the effect might be marginal, we 
note that the temperature gradient signifies that smaller 
and less massive cores are susceptible to externally trig- 
gered collapse. 

5 Conclusions 

The modelling results presented above depend on several 
unknown factors, above all the collisional coefficients of 
H 2 D + and the possible local variations in the ISRF and 
the dust properties. Nevertheless, the observed H 2 D + 
spectrum can be reproduced reasonably well by a self- 
consistent physical model. The strength of H2D + is that it 
selectively traces the dense nuclei of cold cores, and can 
therefore be used to study their internal structure. 

In Oph D we seem to be witnessing a quiescent phase 
of core evolution. The modelling results suggest, how- 
ever, that the balance can be easily disturbed and that the 
core will eventually collapse. The possibility of a tempo- 
rary equilibrium configuration is interesting as it has been 
suggested that Oph D was formed by turbulent fragmen- 
tation (Steinacker et al. 2005). 

The usefulness of H2D + in studies of the early stages 



of star formation has been predicted by astrochemists and 
modellers (e.g. Bergin et al. 2002). Emerging high- 
altitude radio telescopes, most recently APEX, are pro- 
viding empirical tests of these ideas. 
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A Approximation for the collisional 
rates of H 2 D + 

In the calculation of the state-to-state rate coefficients, 
we have adopted the approximation used by Oka & Epp 
(2004) for H3 + H2 collisions. This approximation is 
based on the principles of complete randomness and a 
detailed balance of upward and downward transitions be- 
tween any two levels. This model does not consider any 
nuclear spin restrictions. It is known, however, that in 
cold gas the conversion between the para and ortho states 
of H 2 D + occurs predominantly via collisions with 0-H2. 
The following reactive transitions are possible (Gerlich et 
al. 2002; Walmsley et al. 2004): 

k* 

p-H 2 D + +o-H 2 ^ o-H 2 D + +p-H 2 (1) 



o - H 2 D + + o - H 2 ^ p - H 2 D + + p - H 2 . (2) 

k- 



In molecular clouds the ortho/para ratio of H 2 , o/p- 
H 2 , is likely to be non-thermal (i.e., more o-H 2 than ex- 
pected from thermodynamic equilibrium where o/p-H 2 
~ 9e-™ AK/T ) (Gerlich et al. 2002; Walmsley et al. 2004; 
Floweret al. 2006). Consequently, reaction (A. 1) pumps 
H 2 D + from para to ortho states with the help of the inter- 
nal energy of o-H 2 , and this effect becomes more marked 
at low temperatures where the endothermic "backward" 
reaction is inhibited. 

Because of the small fraction of o-H 2 , the total rates of 
the nuclear spin changing reactions are likely to be much 



lower than those derived from the Oka & Epp approxi- 
mation. We have therefore scaled the Oka & Epp coeffi- 
cients to be consistent with the total rates of the nuclear 
spin-changing reactions (A.l) and (A. 2). The assump- 
tions used in the calculation are summarised below. 

1) The total rate coefficient of rotationally inelastic and 
nuclear spin reactive H 2 D + + H 2 collisions, kc, corre- 
sponds to that given by the Oka & Epp (2004) approxi- 
mation at the same temperature. The total number of col- 
lisions per cm 3 and s is fcc«(H 2 )n(H 2 D + ). The coefficient 
kc can be expanded in terms of state-to-state coefficients, 
Cij, as 



kc 



(3) 



Z(H 2 D + ) 

where ; and / refer to the rotational states of H 2 D + , and Z 
is the rotational partition function. The Oka & Epp coeffi- 
cients dj are normalized so that kc approaches in warm 
gas the Langevin rate coefficient for ion-molecule colli- 
sions. 

We note that our definition of the collisional rate, 
&c«(H 2 ), includes collisions with both p-H 2 and o-H 2 , and 
the rates of both o-p and p-o conversions of H 2 D + . 

2) When collisions dominate, the relative populations 
within the ortho and para states approach those in ther- 
mal equilibrium, although o/pH 2 D + can be non-thermal. 
This is equivalent to the assumption that the rotationally 
inelastic collisions, i.e., 0-0 and p-p transitions, are fast 
compared to nuclear spin changing reactions. 

3) Reactions (A.l) and (A. 2) describe nuclear spin- 
changing collisions summed up over all states, and deter- 
mine o/p H 2 D + . The following conditions are obtained for 
the state-to-state coefficients C p ,- >0 / and C ;, p / correspond- 
ing to para-ortho and ortho-para transitions, respectively: 



Li 8 pie '"' ZfCp^of 
Z(p - H 2 D+) 

Y*i8oi e '"^ "YifC i,pf 

Z(o - H 2 D + ) 



= o/pH 2 fe| + k 2 (4) 
= o/pH 2 ^ + k\ . (5) 



The partition functions Z(p - H 2 D + ) and Z(o - H 2 D + ) 
have a common zero-energy level (the ground state of 
para-H 2 D + ). Equations (A.4) and (A.5) constrain the cor- 
responding sums for the p-p and 0-0 coefficients, C piiP / 
and C ;,o/, through assumption 1). 

We assume that o/p H 2 D + can be obtained from the con- 
dition of chemical equilibrium: 



o/pH 2 D + = 



o/p H 2 k\ + k 2 
o/pH 2 fc+ + k: 



(6) 



The coefficient k 2 is very small at low temperatures, and 
Eq. (A. 6) is essentially the same as Eq. (7) in Gerlich 
et al. (2002). By substitution one finds that a thermal 
o/pH 2 (~ 9 e -i70.4K/r) makes also / p H 2 D + thermal (~ 

9 e -S6AK/Ty 
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The collisional coefficients satisfy the principle of de- 
tailed balance, i.e., n,C,y = n/C/, when the relative pop- 
ulations are given by the Boltzmann distribution. The 
normalization described above ensures that this is also 
true for ortho-para transitions even if o/pH 2 D + is non- 
thermal. 

In practice, all state-to-state collisional rate coefficients 
were initially calculated using the Oka & Epp formula. 
The rate coefficients between the ortho and para states 
were then scaled using (A.4) and (A.5). Thereafter rate 
coefficients within ortho and para states were normalized 
so that assumption 1) is satisfied. 

We used the following values for the total rate coeffi- 
cients in reactions (A.l) and (A.2): k\ = 2 10 9 cm 3 s _1 , 
k~ = k\e- UKIT , k\ = jfc+/18, and k 2 = k+ e - 251K/T (Walm- 
sley et al. 2004, the k 2 coefficients are slightly modified 
to account for the spin statistics). The o/p ratio was as- 
sumed to be 10~ 4 , which is believed to be characteristic of 
a dense dark cloud at an advanced chemical stage (Walm- 
sley et al. 2004; Flower et al. 2006). 

In very dense gas collisional transitions between the 
lowest rotational states of both para- and ortho-H2D + can 
compete with radiative transitions. For example, the ob- 
served ground-state ortho line ( 1 io — » In) is thermalized 
in the core centre (see Fig. 2, middle panel). 

The collisional coefficients derived here have a steeper 
temperature dependence than those used by van der Tak 
(2005) which are proportional to the square root of the ki- 
netic temperature. The (downward) collisional coefficient 
for the 1 io — > In transitions has, however, a relatively 
smooth slope between ~ 7 and ~ 9 K, where the aver- 
age value is ~ 2 10~ 10 cm 3 s~'. This number corresponds 
to the best-fit value that van der Tak (2005) obtained by 
scaling the Black et al. (1990) coefficient by a factor of 
~5. 
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